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Abstract 

We use a method based on the lubrication approximation in conjunction with a residual-based mass-continuity iterative solution 
scheme to compute the flow rate and pressure field in distensible converging-diverging tubes for Navier-Stokes fluids. We employ 
an analytical formula derived from a one-dimensional version of the Navier-Stokes equations to describe the underlying flow model 
that provides the residual function. This formula correlates the flow rate to the boundary pressures in straight cylindrical elastic 
tubes with constant-radius. We validate our findings by the convergence toward a final solution with fine discretization as well as 
by comparison to the Poiseuille-type flow in its convergence toward analytic solutions found earlier in rigid converging-diverging 
tubes. We also tested the method on limiting special cases of cylindrical elastic tubes with constant-radius where the numerical 
solutions converged to the expected analytical solutions. The distensible model has also been endorsed by its convergence toward 
the rigid Poiseuille-type model with increasing the tube wall stiffness. Lubrication-based one-dimensional finite element method 
was also used for verification. In this investigation five converging-diverging geometries are used for demonstration, validation and 
as prototypes for modeling converging-diverging geometries in general. 

Keywords: ID flow, Navier-Stokes, distensible tubes, converging-diverging tubes, irregular conduits, non-linear systems. 


1. Introduction 

The flow of fluids in converging-diverging tubes has many scientific, technological and medical applications such 
as transportation in porous media, filtration processes, polymer processing, and pathological stenoses and aneurysms 
[1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13]. There are many studies about the flow in converging-diverging rigid conduits 
[14, 15, 16, 17, 18, 19, 20, 21] and distensible conduits with fixed cross sections [22, 23, 24, 25, 26, 27, 28] separately 
as well as many other different geometries and fluid and conduit mechanical properties [29, 30, 31]. There is also a 
considerable number of studies on the flow in converging-diverging distensible conduits; although large part of which 
is related to medical applications such as stenosis modeling [32, 33, 34, 35, 36], 

Several methods have been used in the past for investigating and modeling the flow in distensible converging- 
diverging geometries; the majority are based on the numerical discretization methods such as finite element and 
spectral methods although other approaches such as stochastic techniques have also been employed. However, due 
to the huge difficulties associating this subject which combines tube wall deformability with convergence-divergence 
non-linearities, most of these studies are based on substantial approximations and modeling compromises. Moreover, 
they are usually based on very complex mathematical and computational infrastructures which are not only difficult to 
implement and use but also difficult to verify and validate. Also, some of these methods, such as stochastic techniques, 
are computationally demanding and hence they may be prohibitive in some cases. Therefore, simple, reliable and 
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computationally low cost techniques are highly desirable where analytical solutions are not available due to excessive 
difficulties or even impossibility of obtaining such solutions which is the case in most circumstances. 

In this paper we propose the use of the lubrication approximation with a residual-based non-linear solution scheme 
in association with an analytical expression for the flow of Navier-Stokes fluids in straight cylindrical elastic tubes with 
fixed radius to obtain the flow rate and pressure field in a number of cylindrically-symmetric converging-diverging 
geometries with elastic wall mechanical properties. The proposed method combines simplicity, robustness and ease 
of implementation. Moreover, it produces solutions which are very close to any targeted analytical solutions as the 
convergence behavior in the investigated special cases reveals. 

Although the proposed method is related to a single distensible tube, it can also be extended to a network of 
interconnected distensible tubes with partially or totally converging-diverging conduits by integrating these conduits 
into the network and giving them a special treatment based on the proposed method. This approach, can be utilized for 
example in modeling stenoses and other types of flow conduits with irregular geometries as part of fluid flow networks 
in the hemodynamic and hemorheologic studies and in the filtration investigations. 

The method also has a wider validity domain than what may be thought initially with regard to the deformability 
characteristics. Despite the fact that in this paper we use a single analytical expression correlating the flow rate to 
the boundary pressures for a distensible tube with elastic mechanical properties, the method can be well adapted to 
other types of mechanical characteristics, such as tubes with viscoelastic wall rheology, where different pressure-area 
constitutive relations do apply. In fact there is no need even to have an analytical solution for the underlying flow model 
that provides the basic flow characterization for the discretized elements of the converging-diverging geometries in 
the lubrication approximation. What is actually needed is only a well defined flow relation: analytical, or empirical, 
or even numerical [37] as long as it is viable to find the flow in the discretized elements of the lubrication ensemble 
using such a relation to correlate the flow rate to the boundary pressures. 

There is also no need for the geometry to be of a fixed or regular shape as long as a characteristic flow can be 
obtained on the discretized elements, and hence the method can be applied not only to axi-symmetric geometries 
with constant-shape and varying cross sectional area in the flow direction but can also be extended to non-symmetric 
geometries with irregular and varying shape along the flow direction if the flow in the deformable discretized elements 
can be characterized by a well-defined flow relation. The method can as well be applied to non-straight flow conduits 
with and without regular or varying cross sectional shapes such as bending compliant pipes. 

2. Method 

The flow of Navier-Stokes fluids in a cylindrical tube with a cross sectional area A and length L assuming a slip- 
free incompressible laminar axi-symmetric flow with negligible gravitational body forces and fixed velocity profile is 
described by the following one-dimensional system of mass continuity and linear momentum conservation principles 




dA dQ 

t> 0, 

x e [0, L] 

(1) 

dQ 

dt 

8 1 
+ dx ' 

- ° 

t> 0, 

x e [0, L] 

(2) 


In these two equations, Q is the volumetric flow rate, t is the time, x is the axial coordinate along the tube length, 
a is the momentum flux correction factor, p is the fluid mass density, p is the axial pressure which is a function of 
the axial coordinate, and k is the viscosity friction coefficient which is usually given by k = where v is the fluid 
kinematic viscosity defined as the ratio of the fluid dynamic viscosity p to its mass density [38, 39, 40, 27, 41, 42, 21]. 
These relations are usually supported by a constitutive relation that correlates the pressure to the cross sectional area 
in a distensible tube, to close the system in the three variables A, Q and p and hence provide a complete mathematical 
description for the flow in such conduits. 

The correlation between the local pressure and cross sectional area in a compliant tube can be described by many 
mathematical constitutive relations depending on the mechanical characterization of the tube wall and its response to 
pressure such as being elastic or viscoelastic, and linear or non-linear. The following is a commonly used pressure-area 
constitutive elastic relation that describes such a dependency 
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p = f(VI-VT) 

where ft is the tube wall stiffness coefficient which is usually defined by 


(3) 


\fnh 0 E 
1 - q- 


(4) 


A a is the reference cross sectional area corresponding to the reference pressure which in this equation is set to zero 
for convenience without affecting the generality of the results, A is the tube cross sectional area at the actual pressure 
p as opposite to the reference pressure, h a is the tube wall thickness at the reference pressure, while E and g are 
respectively the Young’s elastic modulus and Poisson’s ratio of the tube wall. The essence of Equation 3 is that the 
pressure is proportional to the radius growth with a proportionality stiffness coefficient that is seated by the reference 
area, ft should be remarked that we assume here a constant ambient transmural pressure along the axial direction 
which is set to zero and hence the reference cross sectional area represents unstressed state where A„ is constant along 
the axial direction. 

Based on the pressure-area relation of Equation 3, and using the one-dimensional Navier-Stokes system of Equa- 
tions 1 and 2 for the time-independent flow by dropping the time terms, the following equation correlating the flow 
rate Q to the inlet and outlet boundary areas of an elastic cylindrical tube with a constant unstressed cross sectional 
area over its axial direction can be obtained 


~kL+ Jk 2 L 2 + In (A jn /A ou ) (A]l 2 - A 5 J U 2 ) 

Q = — (5) 

2a In (A in /A ou ) 

where A,-,, and A ou are the tube cross sectional area at the inlet and outlet respectively such that A in > A ou . This 
relation, which in essence correlates the flow rate to the boundary pressures, has been previously [28] derived and 
validated by a finite element scheme. 

The residual-based lubrication approach, which is proposed in the present paper to find the pressure field and 
flow rate in converging-diverging distensible tubes, starts by discretizing the tube in the axial direction into ring-like 
elements. Each one of these elements is approximated as a single tube with a constant radius, which averages the 
inlet and outlet radii of the element, to which Equation 5 applies. A system of non-linear equations based on the mass 
continuity residual and boundary conditions is then formed. 

For a tube discretized into ( N - 1) elements, there are N nodes: two boundaries and ( N - 2) interior nodes. Each 
one of these nodes has a well-defined axial pressure value according to the one-dimensional formulation. Also for the 
interior nodes, and due to the incompressibility of the flow, the total sum of the volumetric flow rate, signed (+/-) 
according to its direction with respect to the node, is zero due to the absence of sources and sinks, and hence (N - 2) 
residual functions which describe the net flow at the interior nodes can be formed. This is associated with two given 
boundary conditions for the inlet and outlet boundary nodes to form N equations. 

A standard method for solving such a system is to use an iterative non-linear simultaneous solution scheme such 
as Newton-Raphson method where an initial guess for the interior nodal pressures is proposed and used in conjunction 
with the Jacobian matrix of the system to find the pressure perturbation vector which is then used to adjust the pressure 
values and repeat this process until a convergence criterion based on the size of the residual norm is reached. The 
process is based on iterative solving of the following equation 


JAp = -r (6) 

where J is the Jacobian matrix, p is the vector of variables which represent the pressure values at the boundary and 
interior nodes, and r is the vector of residuals which, for the interior nodes, is based on the continuity of the volumetric 
flow rate as given by 


fj = = 0 (7) 

/= 1 
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where m is the number of discretized elements connected to node j which is two in this case, and Q, is the signed 
volumetric flow rate in element i as characterized by Equation 5. Equation 6 is then solved in each iteration for Ap 
which is then used to update p. The convergence will be declared when the norm of the residual vector, r, becomes 
within a predefined error tolerance. More details about this solution scheme can be found in [27, 11]. 

3. Implementation and Results 

The proposed residual-based lubrication method was implemented in a computer code and flow solutions were 
obtained for an extensive range of fluid, flow and tube characterizations such as fluid viscosity, flow profile, and tube 
mechanical properties. Five regular converging-diverging axi-symmetric tube geometries were used in the current in- 
vestigation; representative graphic images of these geometries are shown in Figure 1, while the mathematical relations 
that describe the dependency of the tube radius, R, on the tube axial coordinate, x, for these geometries are given in 
Table 1 . A generic converging-diverging tube profile demonstrating the setting of the coordinate system for the R-x 
correlation, as used in Table 1, is shown in Figure 2. These geometries have been used previously [20, 21] to find 
flow relations for Newtonian and power law fluids in rigid tubes. A representative sample of the flow solutions on 
distensible converging-diverging tubes are also given in Figures 3-7. 



(e) Sinusoidal 


Figure 1 . Converging-Diverging tube geometries used in the current investigation. 

In all flow simulations, including the ones shown in Figures 3-7, we used a range of evenly-divided discretization 
meshes to observe the convergence behavior of the solution with respect to mesh refinement. In all cases we noticed 
an obvious trend of convergence with improved meshing toward a final solution that does not tangibly improve with 
further mesh refinement. We also used in these flow simulations a rigid conduit flow model with the same geometry 
and fluid and flow properties where the flow in the rigid discretized elements was modeled by Poiseuille equation. The 
purpose of this use of the rigid model is to assess the solution scheme and test its convergence to the correct solution 
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Figure 2. Generic converging-diverging tube profile demonstrating the coordinate system setting for the correlation between the axial coordinate x 
and the tube radius R used in Table 1 . 


Table 1 . The equations describing the dependency of the tube radius R on the tube axial coordinate x for the five converging-diverging geometries 
used in the current investigation. In all these relations - 1 < x < | and R nun < R max where R m j„ is the tube minimum radius at x = 0 and R max is 
the tube maximum radius at x = ± as demonstrated in Figure 2. 


Geometry 
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(^)_(^)cos( ¥ ) 


because for Poiseuille-type flow with rigid geometries we have analytical solutions, given in Table 2, that correlate the 
flow rate to the pressure drop. Poiseuille-type solutions can also provide a qualitative indicator of the sensibility of 
the distensible solutions; for instance we expect the deviation between the two solutions to decrease with increasing 
the stiffness of the elastic tube. In all cases the correct quantitative values and qualitative trends have been verified. 

Each one of Figures 3-7 shows a sample of the numeric solutions for two sample meshes used for the distensible 
flow geometry alongside the converged Poiseuille-type solution for the given fluid and tube parameters. The reason for 
showing two meshes for the distensible geometry is to demonstrate the convergence behavior with mesh refinement. 
In all cases, virtually identical solutions were obtained with meshes finer than the finest one shown in these figures. 

It should be remarked that in all the distensible flow simulations shown in Figures 3-7 we used a = 4/3 to match 
the rigid Poiseuille-type flow profile [21] which we used, as indicated already, as a test case. However, for the purpose 
of testing and validating the distensible model in general we also used an extensive range of values greater than and 
less than 4/3 for a without observing incorrect convergence or convergence difficulties. In fact using values other 
than a - 4/3 makes the convergence easier in many cases [11], 

An interesting feature that can be seen in Figure 4 is that all the pressure profile curves are almost identical as well 
as the flow rates. The reason is that, due to the high tube stiffness used in this example, the distensible tube solution 
converged to the rigid tube Poiseuille-type solution. A more detailed comparison between the Poiseuille-type rigid 
tube flow and the Navier-Stokes one-dimensional elastic tube flow with high stiffness is shown in Figure 8 where the 
results of Figures 3-7 are reproduced using the same fluid, flow and tube parameters but with high tube stiffness by 
using large /Ts. As seen in Figure 8 the elastic tube flow converges almost identically to the Poiseuille-type rigid 
tube flow with increasing the tube wall stiffness in all cases. This sensible and correct trend can be regarded as 
another verification and validation for the residual-based method and the related computer code. Similar results have 
also been obtained in [41] in comparing the rigid and distensible models for the flow in networks of interconnected 
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Table 2. The equations describing the dependency of the flow rate Q on the pressure drop A p for the rigid tubes with the five converging-diverging 
geometries of Table 1. These relations were previously [21] derived and validated. 
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straight cylindrical tubes. More detailed comparisons between the rigid and distensible one-dimensional flow models 
can be found in the aforementioned reference. 

It should be remarked that the critical value of / 3 at which the distensible flow solution converges to the rigid flow 
solution depends on several factors such as the fluid and flow parameters as well as the geometry of the tube and the 
pressure field regime characterized by the applied boundary conditions at the inlet and outlet where their size and the 
magnitude of their difference play a decisive role. Another remark is that the shape of the pressure profile curve is 
highly dependent on the geometric factors such as j+-, and ratios. It also depends on the fluid and tube 
mechanical properties, such as fluid viscosity and tube wall stiffness, and the magnitude of pressure at the inlet and 
outlet boundaries. 

The opposite to what in Figure 4 can be seen in Figure 5 for the hyperbolic geometry where we used very low 
stiffness and hence the elastic model deviated largely from the rigid model. This also affected the dependency of 
convergence rate on discretization where the discrepancy between the solutions of the coarse and fine meshes was 
more substantial than in the other cases for similar coarse and fine meshes. In general, the deviation between the rigid 
and distensible flow models is maximized by reducing the stiffness, and hence increasing the tube distensibility, while 
other parameters are kept fixed. 

Another interesting feature is that in the flow solution of Figure 6 there is a big difference between the flow rate 
of the elastic and rigid tubes. This can be explained largely by the significant deviation from linearity due to the large 
values of the inlet and outlet boundary pressures, as well as the large size of their difference, with a relatively low 
stiffness. This indicates that the rigid tube flow model is not a suitable approximation for simulating and analyzing the 
flow in distensible tubes and networks, as it has been done for instance in some hemodynamic studies. More detailed 
discussions about this issue can be found in [41], 

In Figure 9 we draw the geometric profile of the elastic tube for the stressed and unstressed states for the five 
examples of Figures 3-7 where we plot the tube radius versus its axial coordinate for these two states. As seen, these 
plots show another sensible qualitative trend in these results and hence provide further endorsement to the residual- 
based method. It is needless to say that in Figures 3-9 the inlet boundary is at x = 0 while the outlet boundary is at the 
other end. 

4. Tests and Validations 

We used several metrics to validate the residual-based method and check our computer code and flow solutions. 
First, we did extensive tests on distensible cylindrical tubes with fixed radius using different fluid, flow and tube 
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Figure 3. Axial pressure as a function of axial coordinate for a converging-diverging elastic tube with conic geometry having L = 0.1 m, R m i„ = 
0.005 m, R max = 0.01 m, and p = 236.3 Pa.m. The fluid properties are: p = 1000 kg.m -3 and p = 0.01 Pa.s while the inlet and outlet pressures 
are: p\ = 1000 Pa and p Q = 0.0 Pa. The Poiseuille-type flow uses a rigid tube with the same unstressed geometry and the same p and boundary 
pressures. The converged flow rate for the elastic Navier-Stokes and rigid Poiseuille-type flows are respectively: Q e = 0.000255889 m 3 .s -1 and 
Q r = 0.000842805 m 3 .s -1 while the analytic flow rate for the rigid tube as obtained from the first equation in Table 2 is Q a = 0.000841498 m 3 .s -1 . 



Figure 4. Axial pressure as a function of axial coordinate for a converging-diverging elastic tube with parabolic geometry having L = 0.013 m, 
R m in = 0.0017 m, Rmax = 0.0025 m, and p = 28059.0 Pa.m. The fluid properties are: p = 1100 kg.m -3 and p = 0.006 Pa.s while the 
inlet and outlet pressures are: p\ = 2000 Pa and p 0 = 1000 Pa. The Poiseuille-type flow uses a rigid tube with the same unstressed geometry 
and the same p and boundary pressures. The converged flow rate for the elastic Navier-Stokes and rigid Poiseuille-type flows are respectively: 
Q e = 6.58209 x 10 -5 m 3 .s -1 and Q, = 6.62929 x 10 -5 m 3 .s -1 while the analytic flow rate for the rigid tube as obtained from the second equation 
in Table 2 is Q a = 6.62051 x 10 -5 m 3 .s -1 . 
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Figure 5. Axial pressure as a function of axial coordinate for a converging-diverging elastic tube with hyperbolic geometry having L = 0.09 m, 
R min = 0.004 m, R m ax = 0.006 m, and [5 = 23.6 Pa.m. The fluid properties are: p = 800 kg.m -3 and p = 0.002 Pa.s while the inlet and 
outlet pressures are: p, = 1500 Pa and p 0 = 500 Pa. The Poiseuille-type flow uses a rigid tube with the same unstressed geometry and the 
same p and boundary pressures. The converged flow rate for the elastic Navier-Stokes and rigid Poiseuille-type flows are respectively: Q e = 
0.000147335 m 3 .s -1 and Q, = 0.000934645 m 3 .s -1 while the analytic flow rate for the rigid tube as obtained from the third equation in Table 2 is 
Q a = 0.000933394 m 3 .s -1 . 



Figure 6. Axial pressure as a function of axial coordinate for a converging-diverging elastic tube with hyperbolic cosine geometry having L = 0.7 m, 
R m i n = 0.05 m, R max = 0.08 m, and j. 3 = 3889.4 Pa.m. The fluid properties are: p = 700 kg.m -3 and p = 0.0075 Pa.s while the inlet and outlet 
pressures are: p t = 2500 Pa and p Q = 700 Pa. The Poiseuille-type flow uses a rigid tube with the same unstressed geometry and the same p and 
boundary pressures. The converged flow rate for the elastic Navier-Stokes and rigid Poiseuille-type flows are respectively: Q e = 0.0427687 m 3 .s -1 
and Q r = 1.4184 m 3 .s -1 while the analytic flow rate for the rigid tube as obtained from the fourth equation in Table 2 is Q a = 1.416296 m 3 .s -1 . 
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Figure 7. Axial pressure as a function of axial coordinate for a converging-diverging elastic tube with sinusoidal geometry having L = 6.5 m, 
R m in = 0.2 m, R ma x = 0.5 m, f) = 5064.2 Pa.m. The fluid properties are: p = 900 kg.m 5 and p = 0.06 Pa.s while the inlet and outlet pressures 
are: p, = 1800 Pa and p„ = 300 Pa. The Poiseuille-type flow uses a rigid tube with the same unstressed geometry and the same // and boundary 
pressures. The converged flow rate for the elastic Navier-Stokes and rigid Poiseuille-type flows are respectively: Q e = 0.396769 m 3 .s~* and 
Q r = 8.74955 m \s 1 while the analytic flow rate for the rigid tube as obtained from the fifth equation in Table 2 is Q a = 8.73370 m\s o 


parameters where the method produced results identical to the analytical flow solutions given by Equation 5. Although 
this test is based on a simple limiting case and hence it may be regarded as trivial, it provides sufficient validation 
for the basic approach and the reliability of the code. We also investigated the convergence behavior, outlined in the 
previous section, as a function of discretization; in all cases it was observed that the residual-based method converges 
to a final solution with the use of finer meshes where it eventually stabilizes without tangible change in the solution 
with more mesh refinement. This convergence behavior is a strong qualitative indicator for the accuracy of the method 
and the reliability of the code. As indicated previously, we used evenly-divided regular meshes in all simulations. 

We also used the discretized Poiseuille-type flow in the same converging-diverging geometry but with rigid wall 
mechanical characteristics to validate the solutions, as discussed in the previous section. As seen, we observed in 
all cases the convergence of the Poiseuille-type solutions on using reasonably fine meshes to the analytical solutions 
with errors that are comparable to the machine precision and hence are negligible as they are intrinsic to any machine- 
based numerical method. Since the elastic and rigid models are based on the same mathematical and computational 
infrastructure, the convergence of the rigid flow model to the correct analytical solution can be regarded as an indirect 
endorsement to the elastic model. The convergence of the elastic model solution to the verified rigid model solution 
with increasing tube wall stiffness is another indirect support for the elastic model as it demonstrates its sensible 
behavior. 

As another way of test and validation, we produced a sample of lubrication-based one-dimensional finite element 
solutions which are obtained by discretizing the converging-diverging distensible geometries and applying the pressure 
continuity, rather than the Bernoulli energy conservation principle, as a coupling condition at the nodal interfaces [27, 
13] to match the assumptions of the residual-based method which couples the discretized elements by the continuity 
of pressure condition [11]. The finite element results were very similar to the residual-based results although the 
convergence behavior was generally different. The residual-based method has a better convergence behavior in general 
although this is highly dependent on coding technical issues and implementation specificities and hence cannot be 
generalized. 

With regard to the comparison between the residual-based and finite element methods, they have very similar 
theoretical infrastructure as they are both based on the same formulation of the one-dimensional Navier-Stokes flow. 
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Figure 8. Comparing the converged Poiseuille-type rigid tube flow (solid) to the converged elastic tube flow with high wall stiffness of the given 
ft (dashed) for the five examples of Figures 3-7. In all the five sub-figures, the vertical axis represents the axial pressure in pascals while the 
horizontal axis represents the tube axial coordinate in meters. The converged numeric flow rate in each case for the rigid and elastic models is 
virtually identical to the corresponding Poiseuille-type analytic flow rate given in Figures 3-7. 


In fact the residual-based method is a modified version of the previously proposed [11] pore-scale network modeling 
method for the flow of Navier-Stokes fluids in networks of interconnected distensible tubes by extending the concept 
of a network to serially-connected tubes with varying radii which represent the discretized elements of the converging- 
diverging tubes. Hence the agreement between the residual-based and finite element methods may not be regarded as 
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(a) Conic 



(c) Hyperbolic 



(b) Parabolic 



(d) Hyperbolic Cosine 



Figure 9. Comparing the elastic tube unstressed radius (solid) to the stressed radius (dashed) as a function of the tube axial coordinate for the five 
examples of Figures 3-7. In all the five sub-figures, the vertical axis represents the tube radius in meters and the horizontal axis represents the tube 
axial coordinate in meters as well. 


an entirely independent validation method and that is why we did not do detailed validation by the lubrication-based 
one-dimensional finite element. 
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5. Comparisons 

As indicated previously, the advantages of the residual-based method in comparison to other methods include 
simplicity, ease of implementation, low computational costs, and reliability of solutions which are comparable in their 
accuracy to any intended analytical solutions based on the given assumptions, as the investigated limiting cases like 
rigid and fixed-radius tubes have revealed. These advantages also apply for the residual-based method in compari- 
son to the lubrication-based one -dimensional finite element method plus a better overall convergence behavior. The 
biggest advantage of the finite element method, however, is its applicability to the transient time -dependent flow and 
more suitability for probing other flow-related one-dimensional transport phenomena such as the reflection and prop- 
agation of pressure waves. Therefore, the lubrication-based one-dimensional finite element could be the method of 
choice for investigating transient flow and wave propagation in distensible geometries until proper modifications are 
introduced on the residual-based method to extend it to these modalities. More details about the comparison between 
the residual-based and finite element methods can be found in [11]. 

The residual-based method, as indicated earlier, can also be used for irregular flow conduits in general with cross 
sections that vary in size and shape and even without converging-diverging feature and regardless of being cylindrically 
axi-symmetric as long as an analytical, or empirical, or even numerical [37] relation between the boundary pressures 
and flow rate on a straight geometry with a similar cross sectional shape does exist. Therefore it can be safely claimed 
that the residual-based method has a wider applicability range than many other methods whose explicit or implicit 
underlying assumptions apply only to restricted types of conduit geometry. 

With regard to convergence, each numerical method has its own characteristic convergence behavior which de- 
pends on many factors such as the utilized numerical solvers and their underlying mathematical and computational 
theory, the nature of the physical problem, the employed convergence support techniques, coding technicalities, and 
so on. Hence it is not easy to make a definite comparison for the convergence behavior between different numerical 
methods. However, we can say that the residual-based method has in general a better rate and speed of convergence 
in comparison to other commonly-used numerical methods. More details about convergence issues and convergence 
enhancement techniques can be found in [11]. 

On the other hand, the residual-based method has a number of limitations based on its underlying physical assump- 
tions, as stated in section 2, as well as limitations rooted in its one -dimensional nature that restricts its applicability 
to modeling axially-dependent flow phenomena and hence excludes phenomena related to other types of dependency. 
However, most of these limitations are shared by other comparable methods. 

6. Conclusions 

A simple and reliable method based on the lubrication approximation in conjunction with a non-linear simultane- 
ous solution scheme based on the continuity of pressure and volumetric flow rate with an analytical solution correlating 
the flow rate to the boundary pressures in straight cylindrical elastic tubes with constant radius is used in this paper to 
find the flow rate and pressure field in distensible tubes with converging-diverging shapes. Five converging-diverging 
axi-symmetric geometries were used for demonstrating the applicability of the method and assessing its merit. 

The method is validated by its convergence behavior with finer discretization as well as comparing the equiv- 
alent Poiseuille-based flow to the analytical solutions which were obtained and validated previously. A sample of 
lubrication-based one-dimensional finite element solutions have also been obtained and compared to the residual- 
based solutions; these results show very good agreement. The method was also tested on limiting cases of elastic 
cylindrical tubes with fixed radius, where it produced results identical to the analytical solutions, as well as the con- 
vergence to the established rigid tube flow with increasing tube wall stiffness. 

The method can be extended to geometries other than cylindrically axi-symmetric converging-diverging shapes as 
long as a flow characterization relation can be provided for the discretized elements; whether analytical or empirical or 
even numerical. The method can also be extended beyond the use in computing the flow in single tubes to computing 
the flow in networks of interconnected distensible conduits which are, totally or partially, characterized by having 
converging-diverging geometries, or variable cross sectional shapes or curving structure in the flow direction to be 
more general. 

Many industrial and medical applications, such as material processing and stenosis modeling, can benefit from this 
approach which is easy to implement and integrate with other flow modeling techniques. Moreover, it produces highly 
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accurate solutions with low computational costs. An initial investigation indicates that its convergence behavior is 
generally superior to that of the traditional numerical techniques such as the one-dimensional finite element especially 
with the use of convergence enhancement techniques. 
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8. Nomenclature 

a correction factor for axial momentum flux 
P stiffness coefficient in the pressure-area relation 
k viscosity friction coefficient 

p fluid dynamic viscosity 

v fluid kinematic viscosity 

p fluid mass density 

Poisson’s ratio of tube wall 
A tube cross sectional area at actual pressure 

Aj n tube cross sectional area at inlet 

A„ tube cross sectional area at reference pressure 

A ou tube cross sectional area at outlet 

E Young’s elastic modulus of the tube wall 
/ flow continuity residual function 

h a tube wall thickness at reference pressure 

J Jacobian matrix 

L tube length 

N number of discretized tube nodes 

p pressure 

p pressure vector 

Pi inlet pressure 

p„ outlet pressure 

A p pressure drop 

Ap pressure perturbation vector 

Q volumetric flow rate 

Q a analytic flow rate for rigid tube 

Q e numeric flow rate for elastic tube 

Q r numeric flow rate for rigid tube 

r residual vector 

R tube radius 

R max maximum unstressed tube radius 

R m j n minimum unstressed tube radius 

t time 
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x tube axial coordinate 
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